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ABSTRACT 

We present results of two- and three-dimensional Particle-In-Cell simulations 
of magnetic-turbulence production by isotropic cosmic-ray ions drifting upstream 
of supernova remnant shocks. The studies aim at testing recent predictions of a 
strong amplification of short-wavelength magnetic field and at studying the sub- 
sequent evolution of the magnetic turbulence and its backreaction on cosmic ray 
trajectories. We observe that an oblique filamentary mode grows more rapidly 
than the non-resonant parallel modes found in analytical theory, and the growth 
rate of the field perturbations is much slower than is estimated for the paral- 
lel plane-wave mode, possibly because in our simulations we cannot maintain 
^ fij, the ion gyrofrequency, to the degree required for the plane- wave mode 
to emerge. The evolved oblique filamentary mode was also observed in MHD 
simulations to dominate in the non-linear phase, when the structures are already 
isotropic. We thus confirm the generation of the turbulent magnetic field due to 
the drift of cosmic-ray ions in the upstream plasma, but as our main result find 
that the amplitude of the turbulence saturates at about SB/B ~ 1. The backre- 
action of the magnetic turbulence on the particles leads to an alignment of the 
bulk-flow velocities of the cosmic rays and the background medium. This is an 
essential characteristic of cosmic-ray modified shocks: the upstream flow speed is 
continuously changed by the cosmic rays. The deceleration of the cosmic-ray drift 
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and the simultaneous bulk acceleration of the background plasma account for the 
saturation of the instability at moderate amplitudes of the magnetic field. Pre- 
viously published MHD simulations have assumed a constant cosmic-ray current 
and no energy or momentum flux in the cosmic rays, which excludes a backre- 
action of the generated magnetic field on cosmic rays, and thus the saturation 
of the field amplitude is artificially suppressed. This may explain the continued 
growth of the magnetic field in the MHD simulations. A strong magnetic-field 
amplification to amplitudes 6B ^ Bq has not been demonstrated yet. 

Subject headings: acceleration of particles, cosmic rays, methods: numerical, 
shock waves, supernova remnants, turbulence 



1. INTRODUCTION 



The prime candidate for the sources of Galactic cosmic rays are shell-type supernova 
remnants (SNRs), at whose forward shocks a Fermi-type acceleration process could accel- 
erate particles to PeV-scale energy, although no firm evidence has been established to date 
for hadron production in SNRs. The acceleration arises from pitch-angle scattering in the 
plasm a flows tha t have systematically different velocities upstream and downstream of the 
shock (jBelllll978l ). More detailed studies show that the acceleration efficacy and the resulting 
spectra depend on the orientation angle of the large-scale magnetic field and on the amphtude 



and characteris t ics of magnetic turbulence near the shock (e.g. iGiacalone fc Jokipiil Il996 
Giacalond l2005l : iPelletier et al.l l2006l : iMarcowith et al.l |2006|). Particle confinement n ear the 
shock is supported by self-generated magnetic turbulence (IMalkov fc Diamondll200ll ). which 
is likely generated in the upstream region. As the flow convects the turbulence toward and 
through the shock, the turbulent magnetic field structure at and behind the shock is shaped 
by the plasma interactions ahead of it. Therefore, detailed knowledge of the properties of 
the turbulence is crucial to ascertain all aspects of the acceleration processes: transport 
properties of cosmic rays, the shock structure, thermal particle injection and heating pro- 
cesses. Of particular interest is the question of how efficiently and with what properties 
electromagnetic turbulence is produced by energetic particles at some distance upstream of 
the forward shocks. If the cosmic rays would drive a tur bulent magnetic field to an ampli 



tude much larger than the homogeneous interstellar field (ILucek fc Belli l2000l : iBell fc Lucek 



200X1 ). part icle acceleration cou 



estimated (ILagage &: Cesarskv 



d be faster and extend to higher energies than conventionally 
IQSsh. although it is unclear what energies could be reached 



fIVladimirov et al.ll2006l : lElhson fc Vladimirovl 120081 ) 



We investigate the properties of magnetic turbulence upstream of the shocks of young 
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SNRs. Cosmic rays accelerated at these shocks form a nearly isotropic population of rela- 
tivistic particles that drift with the shock velocity relative to the upstream plasma. Even 
though the shock may accelerate particles over a wide range of energies, the highest-energy 
particles propagate the farthest from the shock, so at some distance upstream one finds 
only particles of th e highest en ergy, which are predominantly ions. Such a system has been 
recently studied by iBelll ( l2004j . hereafter denoted as B04) with a quasi-linear MHD (Magne- 
tohydrodynamics) approach. Bell noted that, rather than resonant Alfven waves, the current 
carried by drifting cosmic rays should efficiently excite non-resonant, nearly purely growing 
(3? a; = cur 0) modes of mag n etic t urbul ence. The ana l ytical results were reproduced 



elletier et al.l J2006h. who noted the need to 



Zirakashvili et al.l (l2008l. hereafter denoted as 



in studies by iBykov &: ToptyginI (120051 ) and 
carefully account for t he return curren t , and 

Z08). Calculations by iBlasi fc Amatol (120071 ) and iReville et al.l (l2007f ) using the linearized 
Vlasov equation confirmed the predictions of the MHD approximation. iBelll (120051 . hereafter 
denoted as B05) noted that a filamentation of the background plasma must be expected as 
soon as a small current imbalance occurs, which may eventually give rise to a filamenta- 
tion of the cosmic-ray current. MHD simulations described by B04, B05 and Z08 indicate 
a strong magnetic-field amplification following a plasma filamentation that turns approxi- 
mately isotropic in the non-linear stage. The results presented in those papers do not indicate 
that a parallel plane-wave mode is initially observed. 

Here we use Particle-In-Cell (PIC) simulations with a view to explore the properties 
of the magnetic turbulence produced: its geometry, non-linear evolution, saturation, and 
associated particle heating. Our kinetic approach also allows us to address the backreaction 
of the magnetic turbulence on cosmic rays. 

We describe the simulation setup in §2, where we also present a rationale for selecting 
sets of simulations according to the theoretical analyses. In §3 the simulation results are 
presented starting with the main three-dimensional experiment, including a discussion of 
scaling and parameter dependences, and the evolution of the particle spectra. We conclude 
with a summary and discussion in §4. 



2. SIMULATION SETUP 

2.1. Simulation Method and Initial Conditions 

Our discussion of the production of magnetic turbulence by the non-resonant streaming 
instability is based on three-dimensional PIC experiments. However, performing large-scale 
simulations in three dimensions poses serious computational demands. Simulations of two- 
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dimensional systems can address problems that involve a wider range of spatial scales, or 
growth on timescales long compared with the electron plasma periods, and therefore our 
three-dimensional experiments are complemented by a series of two-dimensional simula- 
tions. The two-dimensional runs allow us to study the evolution of magnetic turbulence 
at timescales for which the turbulent field structures become larger than the extent of the 
three-dimensional simulation box, and so they enable us to ascertain the influence of the 
limited size of the simulation box on the results of the 3-D experiments. Furthermore, the 
two-dimensional runs permit us to study how these results depend on the choice of physical 
parameters of the simulations. 



The code used in this study is a modified parallel version of TRISTAN (lBunemanlll993l ) . 
a three-dimensional fully relativistic electromagnetic Particle-In-Cell code in Cartesian co- 
ordinates. The modifications to the code include the use of a new fi r st-ord er algorithm for 



the charge-conserving current deposition developed by lUmeda et al.l (120031 ). and the imple- 
mentation of a new method of digital filtering of electric currents. The filtering method uses 
the set of 27-point-average isotropic smoothing and compensation filters, which attenuate 
the short-wavelength noise, increase the overall accu racy, and at the same tim e improve 



agreement with theory at long wavelengths (see, e.g., iBirdsall fc Langdoru l2005l ) . The two- 
dimensional simulations have been performed using the same three-dimensional code, in 
which the grid size in one dimension was restricted to contain three cells only. We thus 
follow the temporal evolution of all three vector components of particle velocities and the 
electromagnetic fields, and can therefore accurately simulate circularly polarized waves as 
found in the analytical theory, requiring only that the wavevector lies in the simulation plane. 
However, a filamentation mode may appear somewhat different in two dimensions. 

In the simulations, an isotropic population of relativistic, monoenergetic cosmic-ray ions 
with Lorentz factor '-jcr = 2 and number density Ncr drifts with Vsh = 0.3c relative to the 
upstream electron-ion plasma and along a homogeneous magnetic field B\\q, carrying a current 
density jcn = eNcRVsh- The ions of the upstream medium have a thermal distribution with 
number density Ni, in thermal equilibrium with the electrons. The electron population with 
thermal velocity Vg^th = 0.002c and density = + Ncr contains the excess electrons 
required to provide charge-neutrality and also drifts with Vd = VghNcR/Ne with respect to 
the upstream ions, so it provides a return current jret = —^NeVd to balance the current in 
cosmic rays. We assume that the return current is entirely carried by the electron component 
of the background plasma because that will most quickly react to an electric field induced 
by charge separation due to the drifting cosmic-ray component. The simulations have been 
performed using computational grids with periodic boundary conditions (see Table 1), which 
keeps the particle number constant throughout the simulations at a value corresponding to 
an average total density of 16 particles per cell. A two-dimensional test simulation with 160 
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particles per cell has also been run to verify that the number of particles used in our studies 
does not significantly influence the results. 

Cosmic-ray ions drift in the — x-direction, antiparallel to the homogeneous magnetic 
field B\\Q. For the two-dimensional simulations we render the 2;-coordinate ignorable, so the 
direction of the drift is in the simulation plane. The electron skindepth Xse — c/ujpe — 
7A, where Upe — {N^e^ /meeoY^"^ is the electron plasma frequency and A is the grid cell 
size. The parameter combination is chosen both to resemble physical conditions in young 
SNRs and, according to the quasi-linear calculations by B04, to be favorable for the rapid 
excitation of purely growing, short-wavclcngth (as compared with cosmic-ray ion gyroradii 
TcRg) wave modes, although we typically cannot maintain the condition uj <^ Vt^. For a 
cold ambient plasma and for wavevectors k\\ parallel to Syo, the dispersion relation in the 
unstable wavevector regime 1 < k\\rcRg < C'^sh/'^l reads 

- v\kl + Cvln = 0, (1) 

rcRg 

for which purely growing solutions arc found. The maximum growth rate '^max at the 
wavenumber k\\max can be determined as 

2 v\ rcRg 2 va rcRg 



Here, rcRg — mjC-\/7^^ — l/eBp is the gyroradius of the (lowest-energy) cosmic rays, va — 
[B'^^q/ lj,Q{Neme + A^jmj)]"^/^ is the Alfven velocity, and ( = Ncrc^^^cr ~ '^/^i'^sh determines 
the strength of the cosmic-ray driving term. From Eq. (1) one can see that the non-resonant 
modes are strongly driven when C'^1h/v\ ^ 1. For the given parameters, this condition 
depends on the magnetic field strength and the ion-electron mass ratio, which also determine 
the wavelength \max = '^n/kmax of the most unstable mode. This wavelength must be much 
smaller than the size of the computational box to enable studies of the evolution of the 
turbulence towards larger scales, and at the same time it must be clearly separated from the 
physical scales in the background plasma. Restricted by these computational requirements, 
for the main three-dimensional experiment (run A) we assume a density ratio NcR/Ni=l/3 
and a reduced ion-electron mass ratio mi/rrie = 10, which gives the ion skindepth Xgi = 
25. 3A and sets the wavelength of the most unstable mode to Xmax = 50A. With this 
choice we have further specified Cv'jf^/v\ ~ 715, the Alfvcnic Mach number Ma = Vgh/vA = 
19.1, and the ratio ujpe/fle = 22.1, where the electron cyclotron frequency fl^ = eB/rrie, 
which corresponds to weakly magnetized conditions of the interstellar medium. The periodic 
boundary conditions allow us to follow the gyration motion of cosmic rays, even though the 
cosmic-ray ion gyroradii, rcRg — 2672A, are much larger than the size of the simulation box, 
which is 19.8x6.1x6.1 in units of X^n-r- 
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2.2. Simulation Pcirameters 

Table 1 compares the parameters and main results of all runs in two and three dimen- 
sions with those of the initial three-dimensional experiment. The follow-up simulations were 
mainly used to probe systems with more realistic mass ratios up to mi/rrie = 500, better sep- 
arated spatial scales Xgi and Xmax, smaller density ratios Ncn/Ni, and also to allow the box 
size to be a higher multiple of the wavelength of the most unstable mode. Run A* uses the 
parameters of run A in a large two-dimensional box of size 50x18 Xmax and is performed to 
confirm that the results of two-dimensional runs are similar to those of the three-dimensional 
simulations. Runs B to E use larger ion-electron mass ratios, whereas runs G to J are for 
smaller density ratios. 

Note that in runs B to E the ion skindepth Xgi increases with mass ratio as {rrii/meY^'^ , 
and thus Xmax needs to be readjusted to separate the plasma and turbulence scales. The 
wavelength of the most unstable mode is thus increased by increasing the strength of the 
homogeneous magnetic field B\\q by approximately the same factor {mi/m^.Y^'^ . In this way 
the ratio Xmax/Xsi and the Alfvcn Mach number Vsh/vA of the cosmic-ray drift of runs A 
and A* are preserved in those runs, and also the dynamic range in wavelength, as given by 
CVgfi/v\, remains unchanged. However, the maximum growth rate 'jmax scales inversely with 
{mi/meY^^, and therefore the three-dimensional studies would clearly require an enormous 
computational effort, if more realistic mass ratios were assumed. We restrict the duration 
of these simulations (runs B and C) to cover only the initial stage of the instabilities, and 
the full non-linear evolution is investigated by means of two-dimensional simulations (runs D 
and E). Note that when the ion-electron mass ratio increases, the upstream plasma becomes 
more magnetized (the frequency ratio Upe/^le decreases; see Table 1), and the cosmic-ray 
gyroradii grow. For run E, in which mi/rrie = 500, these parameters are ujpe/flg — 3.3 and 
rcRg — 19722A. Finally, in run F the wavelength of the most unstable mode has been 
increased to X^ax — 150A to allow for a better separation from the plasma scales. This is 
done by increasing the strength of the regular magnetic field B\\q, whereas the ion-electron 
mass ratio is kept the same as in run A. This means that in this case the Alfven Mach 
number decreases to Ma = 6.3 and also the width of the unstable wavevector range, here 
Cv^fjv\ ~ 80, is an order of magnitude smaller compared with the other runs. For the 
parameters of run F, the ratio cUpe/fle = 7.4 and rcRg — 890A. 

The runs G to J have been performed for smaller values of the density ratio Ncr/N^ 
down to 1/30, which implies a smaller drift velocity of the electrons in the background 
plasma. This reduces the impact of an initial electrostatic, Buneman-type instability, so 
less heating will occur initially and the plasma temperature should be lower throughout the 
simulation. The mass ratio is mi/rrie = 50 in these simulations {Xgi = 51. 9A), which cover a 
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range in plasma magnetizations with parameters ranging from M4 = 13 and Xmax = 500A 
to Ma = 65 and Xmax = 100 A. In all simulations the time step 6t = 0.0715/ct;pe- 



3. RESULTS 

3.1. Three-dimensional Experiment 

The temporal evolution of the magnetic and electric fields and the particle kinetic- 
energy densities in the three-dimensional simulation (run A) is shown in Fig. 1. An initial 
(up to t ~ Imax) f^st growth of the turbulent fields is caused by a Buneman-type beam 
instability between drifting electrons and the ambient ions, which leads to plasma heating 
until the electron thermal velocities become comparable to the electron drift speed Vd (recall 
that the electrons in the background plasma must slowly drift to initially compensate the 
cosmic-ray current), found in the analytical theory In real SNRs the density ratio Ne/NcR 
is likely much larger than the sonic Mach number, and therefore the electron drift speed is 
smaller than their thermal velocity. The initial instability is thus solely a consequence of the 
initial conditions in our simulations, in particular the initial temperature of the background 
medium and the density contrast between cosmic rays and the background plasma. 

At t ~ 2.5 7~^2. wave modes start to emerge that are excited by the cosmic-ray ions 
streaming in the upstream plasma, and turbulent magnetic field is seen mainly in the compo- 
nents transverse to the cosmic-ray drift direction. A plane-wave Fourier analysis is shown in 
Fig. [2] for one component of the magnetic field and in Fig. [3] for the density of the background 
electrons. All Fourier power spectra are averages over the simulation box to reduce the noise. 
Note that for a plane wave of wavelength A propagating at an angle 9 to the drift (— x) di- 
rection, the one-dimensional Fourier spectra will reveal a signal at Ay = A^ = A/ cos^ and 
at A_L = A2 = A/ sin^, respectively. Assuming a plane wave as in the analytical calculations 
of B04 we therefore conclude that the dominant wave mode is oblique at an angle of about 
6 ^ 80° to the x-direction, and has the perpendicular wavelength component A_l ~ 50A, 
which is numerically similar to Xmax- The character of the mode is thus different than pre- 
dicted by the quasi-linear analysis (B04, BOS), which indicated that initially the most rapid 
growth should occur for wavenumbers parallel to B\\o (see also Eq. 1). Here we observe a 
modified filamentation instability, which appears to be faster than the non-resonant par allel 



modes and was also observed in MHD simulations (B04, BOS, ZOS lReville et al.l (j2008af )) to 
dominate in the non-linear phase. At that time, when 6B > B\\q is already established, the 
field structures in the MHD studies have the same spatial scales parallel and perpendicular 
to the cosmic-ray drift direction, as evident from Figure 3 in Z08. As we describe in more 
detail below, the same is true for the magnetic structures in our simulations, which at that 
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point more resemble cavities and isolated peaks. The appearance of an oblique filamen- 
tary mode is in line with earlier results for electron beams which show that the strongest 
growth is typical l y observed not for the p erpendicular case {k ± Vsh), but at a finite angle 



( Bret et al. 2005 : Dieckmann et al. 2008h. beca use there is a cumulative effect of instabil- 



ities that operate in parallel (ILazar et al.l 120071 ). which naturally is not captured in linear 
analytical analysis. The appearance of magnetic power spectra is then further complicated 
by the fact that filamentation is not exactly a transverse mode in the sense that k ■ E ^ 
( iBret et al.ll2007l ). As seen in Fig. [2], the scale in of the dominant feature in the spectrum 
grows as the turbulence develops, but the parallel scales (Ay = A^; ~ 250A) remain roughly 
constant, so after about 207~^2. the structures have similar size parallel and perpendicular 
to the drift. A comparison with Fig. [1] reveals that at those late times 6B± is of the same 
order as -Byo, so the system becomes magnetically nearly isotropic. 

The spatial structure of the turbulent magnetic field and the density of the ambient 
plasma are shown in Figs. HI [5l and[6]for t ^ 10 'Jmaxy found in the analytical theory Fig. [7] 
presents snapshots of the temporal evolution of these structures in the plane perpendicular 
to the cosmic-ray drift direction. To be noted from Fig. H] is that the density variations of 
the ambient electrons and ions are nearly cospatial, apart from some variations on small 
scales that arise largely from statistical fluctuations, so it is sufficient to only show the ion 
density in Figs. [5], O and [71 The background plasma forms filamentary structures of plasma 
voids surrounded by regions of enhanced plasma density. At any given time, the thickness of 
these structures, or their separation, corresponds to the A^ component of the dominant wave 
mode, and the correlation length along the cosmic-ray streaming direction is equal to Ay. 
The cosmic-ray distribution remains homogeneous throughout the duration of the numerical 
experiment, but the voids become completely depleted of ambient ions and electrons except 
for the excess electrons which neutralize the charge in cosmic-ray ions. The currents carried 
by cosmic rays are therefore no longer neutralized, and a strong net current fiows inside 
the filamentary ambient-plasma voids, resulting in magnetic-field lines circling around the 
center of the cavities according to Ampere's law. Such filamented current structures with 
azimuthal magnetic fie lds and radial electric fields are also observed in other simulations 
dNishikawa et al.lbooei ). 



The perpendicular magnetic field is thus concentrated around regions of low background 
plasma density, but constant cosmic-ray density. In the MHD description, the forma- 
tion of cavities in the plasma is due to the jret x B force, which accelerates the ambient 
plasma away from the cente r of the voids, thus causing the cavities to expand (B04, BOS, 
Milosavljevic fc Nakarl l2006l ). In our kinetic simulations, mainly the electrons are subject 
to this force, since they drift and thus contribute to jret- Ambient-ion motion away from 
the cavities, and a restoring infiuence on the electrons, is caused by an electrostatic charge- 
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separation field, which is stronger than the original jret x B force for non-relativistic drifts. 



3.2. Theoretical Analysis of Current Filament Structure 

The formation of the low-density filaments in the ambient plasma can be understood 
with a simple toy model. We use cylindrical geometry with azimuthal symmetry and assume 
the cosmic rays to drift homogeneously along the 2;-axis, thus providing a current density 
JcR = jcR,z = eNcRVsh- As in our simulations, the ions in the background plasma are at 
rest and homogeneously distributed with density Ni, whereas the background electrons are 
initially homogeneously distributed with density A^^oe and drift along the z-axis with velocity 
■^d = "^d.z = VshNc^/NQe to cauccl the cosmic-ray current and charge. Now suppose a small 
displacement of electrons occurs at some location, so that 



5{r — ri) 3- 5{r — 

r2 



(3) 



where we assume for simplicity r2 > ri and introduce r/ > as a small perturbation pa- 
rameter. The displacement breaks the current balance, so that an excess current exists 
with 



6{r — ri) 6{r — 

r2 



(4) 



This current creates an azimuthal magnetic field that is determined by Ampere's law. The 
field is non-vanishing only between ri and r2 and has the amplitude 

B^{r) = fiQ jcR r]— ri <r <r2 (5) 
r 

The small displacement of electrons also implies a violation of charge balance, leading to a 
radial electric field of amplitude 

^2 

Er(r) = ri — ri < r < r2 (6) 

The electric force on the electrons is therefore a factor /v\ ^ 1 larger than the average 
magnetic force, so they are accelerated inward. The ambient ions see on average only the 
electric force, because at least initially they do not drift, which accelerates them outward, 
albeit with an acceleration much smaller than that of the electrons. The cosmic rays are 
virtually unaffected on account of their large random velocities, which implies that they see 
any acceleration only for a very short time. The dominant electric field prevents a separation 
of ambient ions and electrons, and in fact we see in Fig. H] that these two particle species are 
nearly cospatial. The ambient ions and electrons readjust their spatial distributions until 
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charge balance is re-established, but as a bulk they move slightly outward, so the currents 
are not balanced. Consequently, a current-density profile is established that resembles that 
described in Eq. |H although the parameters t], ri, and r2 will not be the same as before. 
Nevertheless, a magnetic field similar to that in Eq. [S] results, that is not compensated by 
an electric field. The v x B force will pull electrons outward and with them the ambient 
ions, because they are electrostatically bound to the electrons. Reusing Eq. H] with modified 
parameters t]', r[ and with Upi as the ion plasma frequency, the radial acceleration of the 
ambient plasma can be written as 



SO the growth time can be estimated as 



ament — 




T]' Wpi NcR Vsh 



which recovers the scaling of Eq. 2 and is identical to the result of BOS for the filamentary 
mode. It also shows the same scaling as the expected growth rate of the parallel plane-wave 
mode (see Eq. [2]). For efficiently accelerating SNRs tmamont could be of the order of hours or 
days, if i]' were close to unity. 

If the initial displacement of the electrons were inward (r^ < in Eq. [3]), an inward 
transport of ambient plasma would result. However, the electrostatic restoring force (Eq. [6]) 
would be stronger on account of the 1/r-profile, and therefore a weaker current imbalance 
would result after the re-establishment of charge balance in the background plasma, thus 
imposing a preference for outward displacement of the ambient plasma and the creation of 
plasma voids. 



3.3. Non-linear Evolution of Current Filaments 

The expansion of the plasma cavities is visible in Figure U\ and also as the evolution 
toward larger scales in the power spectra shown in Figs. [2] and [31 The growth of the cavities 
and the subsequent merging of the adjacent plasma voids lead to a compression of the plasma 
between cavities and also to an amplification of the magnetic field. Because the magnetic 
field has a preferred orientation around each cavity, the magnetic field lines may cancel each 
other in the space between the voids (Fig. The turbulence is almost entirely magnetic 
and there is no large-scale turbulent electric field structure accompanying the filamentary 
distribution of ambient plasma and perpendicular magnetic field. The electric fields result 
from thermal plasma motions, the level of which can be estimated from Figure HI 
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Compared with Bell's MHD studies, the growth of the magnetic perturbations in our 
kinetic modeling is much slower than calculated for the parallel planar mode. The initial 
growth rate of the perpendicular-field turbulence is only ~ 0.2 7^0^ prior to t'jmax ~ 7, 
and becomes smaller during the later evolution, possibly after nearby plasma cavities have 
started to merge. In the sense of the toy model discussed in §3.2 (Eq. [8]), the perturbation 
parameter must be small, 77' < 0.1, or neighboring filaments partially cancel the magnetic 
field between them. Moreover, the amplitude of the turbulent component never considerably 
exceeds the amplitude of the regular field. Already after about 207~^^, when SB± f» Bp, 
the growth of the turbulence starts to saturate. As shown in Figures [1], [21 and [71 the initial 
filamentary structure of plasma cavities and surrounding magnetic fields becomes disrupted 
at this stage, and the turbulence is nearly isotropic and highly non-linear. In the simulation 
frame, which is the initial rest frame of the ambient ions, the turbulent field structures start 
to move as a consequence of being embedded in the background plasma which itself starts 
to drift (see § 13. 5p . Merging of the magnetic structures leads to further amplification of 
the magnetic field through compression, but at very slow rate. The peak amplitude of the 
average turbulent field, 5B'^°-^ ~ 3.5i?||o, is reached at t-^max ~ 40, after which the field starts 
to dissipate. This is related to the increasing length scales and dissipation of the structures 
in the ambient-plasma density as evident from Figure [3] The electron and ion distributions 
become homogeneous, but still evolve nearly cospatially, so that currents produced in this 
process are weak and thus the gradual decay of the magnetic field is slow. 

Also clearly visible in Figure [3] is that at t > 257~^^ the dominant structures become 
larger than the size of our simulation box. Using two-dimensional simulations performed on 
a large computational grid (run A*), we have verified that the long-time evolution of the 
magnetic turbulence and the plasma density fluctuations in the three-dimensional experiment 
are not considerably influenced by the size of the simulation box. Figure [8] shows the 
temporal evolution of the energy densities in particles and fields for run A*, and Figure [HI 
presents the power spectrum evolution of the perpendicular magnetic field component. To 
be noted from the figures is that the initial growth of the perpendicular magnetic field, 
associated in part with the Buneman-type instability, leads to the higher initial amplitude 
of the transverse field component in the two-dimensional simulation, and thus 5B^ reaches 
the amplitude of the regular field faster compared with run A. Otherwise, the results of the 
three-dimensional experiment are qualitatively very well reproduced in the two-dimensional 
run. At about t'-^max ~ 35 the peak amplitude of the average turbulent field is reached, and 
its value, dB"^"-^ ~ 3.1i?||o, is very close to that obtained in the three-dimensional simulation. 
Also, the wavelengths (A^, Ay) of the dominant wave mode of the magnetic-field turbulence 
are numerically similar to those in run A. The dominant turbulent field structures during 
the late-time evolution are well contained inside the simulation box (Fig. [9]), so that the 
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eventual dissipation of the magnetic field is not affected by the boundary conditions. 



3.4. Scaling and Parameter Dependence 



The question arises to what extent the parameter choice in our simulations influences 
the plasma dynamics, in particular the turbulence growth rates. We have used the zero- 
temperature kinetic calculations of iBlasi fc Amatol (120071 ) to verify that our choice of the 
reduced ion-electron mass ratio mi/rrie = 10 and a monoenergetic cosmic-ray spectrum with 
particle Lorentz factor 7cr = 2 has no impact on the growth rate and spatial scale of the 
instability. For that purpose we have re-evaluated their approximations for u '€i fl*, the 
nonrelativistic ion gyrofrequency, and as the only impact of our parameter choice we find an 
addition al term in the di s persio n relation that scales with the mass ratio. In the notation of 
Eq. 2 in lBlasi fc Amatol (120071 ). the dispersion relation would read 



7 2 ""e 
VaK 

rrii 



Nr 



N,. 



^sh 



(9) 



For a realistic mass ratio and the parameters used by lBlasi fc Amatol (120071 ). the additional 
term is negligible, but for the standard parameters in our simulations it is not. In fact, 
the left-hand side of Eq. [9] is negative for all k, and hence for one polarization we find a 
nearly purely growing solution with a growth rate slightly larger than for irti/me = 1836 
for all krcRg ^ 1 because the functions Ji ^ 1 and /2 ^ 1 for those wavenumbers. These 
calculations suggest that the small growth rate of the non-resonant streaming instability in 
our simulations is not caused by the reduced ion-electron mass ratio, in line with the results 
of additional 2-D simulations using more realistic mass ratios (runs B-E), whose results are 
described below. 

Figure [10] shows the temporal evolution of the energy density in electromagnetic fields 
and particles for a two-dimensional run with rrii/ me = 500 (run E). One can note from the 
figure and also from Table 1 that the fundamental properties of the magnetic turbulence 
observed in the three-dimensional experiment (run A) are also seen in the additional simula- 
tions for mi/rrif. = 40, 100, 150 and 500. In particular, for all mass ratios the dominant wave 
mode, which is the oblique filamentary mode, has approximately the same growth rate, much 
slower than predicted by the quasi-linear estimates, and a wavevector with the same incli- 
nation to the drift direction 6 ^ 70°. The two-dimensional simulations for high ion-electron 
mass ratios (runs D and E) also show that the saturation level of the magnetic turbulence 
is similar to that observed in the three-dimensional experiment with mj/mg = 10. 

Finite-temperat ure effects in background plasma can limit the growth rate of the non- 
resonant instability (iReville et al.ll2006l . l2007l ). Since our simulations do not exactly repro- 
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d uce the analyt i cal est imates for the zero-temperature limit, we used the analytical scalings 



of iReville et al.l (120071 ) to assess the role of the thermal effects. The initial heating of the 
ambient-particle populations associated with, e.g., the Buneman instability leads to ion tem- 
peratures of the order B = /c^r/mjC^ ~ 10"'* — 10~^. At these temperatures the growth rates 
for the most unstable wave modes should be reduced as 7maa;(0 = 10"**) ~ 0.9 'jmaxiQ = 0) 
and 7max(0 = 10~^) ~ O.b'frnaxi'S) = 0)^ This means that thermal effects alone cannot 
explain why in our simulations a filamentation mode with a growth rate of O.2 7max(0 = 0) 
is observed, but not the parallel non-resonant mode. 

For the two-dimensional run F with mi/me = 10 we have increased the strength of the 
regular magnetic field component i?||o. In that case the turbulence scales (Xmax) are much 
better separated from the plasma scales; in fact, Xmax is a factor of three larger than in 
the runs A and A*, whereas the electron and ion skindepths are the same as before. In 
this way we also probe the production of magnetic turbulence in more magnetized upstream 
plasma, in which the shock is weaker with an Alfvenic Mach number Ma = 6.3. The 
results are presented in Figures [11] and [12] and summarized in Table 1, where the reader 
will note that the growth of the turbulence is faster than that for stronger shocks with a 
maximum rate 0.35 'jmax- However, the saturation level is essentially the same as in all the 
other simulations. Furthermore, the faster evolution allows us to better follow the late-time 
behavior of the system, and we clearly observe the turbulent magnetic field to slowly decay 
until at the end of the simulations the energy density in the turbulent magnetic field is within 
a factor 2 of that in the homogeneous magnetic field. We conclude that there is no evidence 
for significant magnetic field amplification, although we do observe a small change in the 
dominant oblique filamentary mode, which has a wavevector that is slightly better aligned 
with the drift direction {9 = Z{k, Vg) ~ 53°), which appears to be a general trend when Xmax 
is expected to be a higher multiple of the ion skindepth. 

We also performed a number of simulations with reduced Ncn/Ni, for which the ex- 
pected growth rate of the non-resonant instability decreases according to the quasi-linear 
estimates. A smaller cosmic-ray current implies that the electrons in the background plasma 
can compensate that current with a smaller drift velocity, thus lessening the impact of the 
initial electrostatic Buneman-type instability and the associated heating. The runs G to 
J therefore involve lower plasma temperatures throughout the duration of each simulation. 
The growth of the magnetic turbulence becomes very slow in these simulations, which use 
the same 6t <^ uo^^ as runs A-F, and therefore in run J with Ncn/Ni = 1/30 we follow only 
the initial evolution. The runs G-I are for Ncr/Ni = 1/10 and three different magnetiza- 



^We are indebted to B. Reville for providing these estimates to us. 
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tions, implying that the expected turbulence scales, Xmax, are between 2 times and 10 times 
the ion skindepth. The mass ratio is mi/rrie = 50 for all runs with small density ratios. As 
in the case of runs A-F, the mean amplitude of the turbulent magnetic field peaks at a few 
times the homogeneous field strength and we observe an oblique filamentation mode, not the 
parallel mode expected from quasi-linear theory. 



3.5. The Evolution of the Particle Spectra 

In all runs a random sample of 0.1% - 0.3% of each particle species (electron, ion, cosmic 
ray), corresponding to about 10^ particles, is selected every few hundred timesteps to observe 
the evolution of the velocity distributions. The average particle velocity is the particles' drift 
velocity and is plotted for all three species in Fig. [13] for run A. To be noted from the figure 
is the disappearance of a relative drift between the cosmic rays and the background ions 
and electrons: after approximately 40 7~^^ all particles drift with approximately the same 
bulk velocity of about 0.14 c. At the same time the turbulent magnetic field reaches its 
peak energy density, as is shown in Fig. [TJ The same basic behavior is seen in the other 
simulations as well. This finding is in line with the notion that the cosmic-ray drift relative 
to the background plasma drives the turbulence: as the relative velocity between cosmic rays 
and plasma decreases, so does the source of the non-resonant streaming instability. We have 
performed a test simulation (run I* in Table 1), in which we do not allow any acceleration 
of cosmic rays by increasing their particle mass by a factor 2 x 10^ and leaving all other 
parameters as in run I. Indeed, the cosmic-ray drift remains constant and the magnetic field 
grows to a bit more than twice the amplitude as compared with run I. The background plasma 
would still experience a bulk acceleration, but now the magnetic-field growth terminates only 
when the plasma has assumed the full original cosmic-ray drift speed rather than about half 
of it as in all the other simulations, thus delaying the saturation. 

The bulk acceleration of the background plasma plays a key role in shaping cosmic-ray 
modified shocks, although in realistic SNRs the process must be expected to operate on a 
timescale similar to that on which the local conditions change on account of the approaching 
shock. In a steady state the upstream plasma will assume a bulk velocity that is essentially 
determined by the local density of cosmic rays. In any case one expects a continuous change 
in the plasma flow velocity, which permits the overall compression ratio between far upstream 
and far downstr eam to be much larger t han the limit given by the Rankine-Hugoniot jump 



conditions (e.g. IVladimirov et al.l l2006l ). Our simulations detail how the changes in the 



bulk-flow properties relate to the average amplitude of the magnetic field. 



Besides the buildup of magnetic turbulence and the bulk acceleration of the background 
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plasma we can also study the heating of the plasma, which is of particular interest for models 
of cosmic-ray acceleration because it limits the sonic Mach number of the shock and hence 
could reduce the acceleration efficiency. After a Lorentz transformation into the bulk-flow 
frame of a given plasma component (see Appendix |X] for details), measurements of the 
azimuth-integrated particle distribution can reveal anisotropies and spectral evolution. Such 
an analysis shows that all particle species, cosmic rays as well as the background ions and 
electrons, remain moderately isotropic in their instantaneous bulk-flow rest frames, but a 
certain stretching of the distributions along the drift direction is clearly seen. Therefore, 
care must be exercised in deriving and interpreting the momentum spectra of the particles. 

The cosmic rays start the simulation with an isotropic and monoenergetic distribution 
with Lorentz factor 7cr = 2. During the simulation isotropy is maintained to within 0.5%, 
and the momentum distribution of the cosmic rays continuously broadens, but remains 
approximately a Gaussian centered on mc^J 7^j^ — 1 with a mode that does not change 
appreciably. 

The momentum distributions of the background ions and electrons intermittently de- 
velop a high-energy tail, at least part of which can be attributed to anisotropy because 
the particles in the high-energy tail predominantly move along the drift direction. A more 
detailed inspection, however, reveals that most of the anisotropies arise from the initial 
electrostatic Buneman-type instability, because they appear long before the non-resonant 
magnetic instability sets in. The initial electrostatic instability involves essentially only elec- 
tric flelds parallel to the drift direction which signiflcantly stretch the distributions of slow 
particles, which is why the background ions are most strongly affected with peak anisotropies 
of 20%. The Lorentz force contributed by the perpendicular components of the magnetic 
flelds provides re-isotropization which is faster for electrons than for ions on account of the 
difference in gyrofrequency, but even in the case of the ions the anisotropy is down to < 3% 
after about t7max — 15 in the main three-dimensional simulation. 

As outlined in appendix \^ we can split the kinetic energy density of particles into the 
components associated with the bulk motion and the random motion. We denote increases 
in the random energy density as heating, but the reader should note that this does not 
necessarily imply a thermalization. Figure [T] shows the temporal evolution of the random 
energy density for all three particle species in run A. Substantial heating of both the electrons 
and the ions in the background plasma is observed early in the simulation as a result of the 
initial electrostatic instability. Further strong heating occurs after )f:7niax — 20, when the 
random energy density of all background plasma species is always higher than their bulk 
energy density. The mean random kinetic energy per particle, or "temperature," is typically 
different for the electron and ions in the background plasma. After the initial Buneman-type 
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instability the electrons are hotter, but roughly at the time the filamentation instability has 
taken over and turned non-linear, the ions are strongly heated and they eventually become 
hotter than the electrons. These results on the anisotropy and heating are seen in the 
main three-dimensional simulation, but qualitatively the same behavior is observed also in 
two-dimensional runs. 

If the observed heating of the background plasma were real, then the upstream medium 
of an SNR forward shock would be so hot that the forward shock itself could be only very 
weak, if it existed at all. There are reasons to assume that our simulations for technical 
reasons overestimate the heating, and therefore we cannot make firm statements on the 
temperature of a real plasma upstream of SNR shocks with efficient particle acceleration. 
One reason is that we must limit the number of simulated particles to simulate the long- 
term evolution in large simulation boxes. Statistical fluctuations then give rise to small-scale 
electric fields that heat the plasma. A short test simulation with large particle number 
(160 per cell in total) was performed and indeed the plasma temperature was observed to 
be somewhat reduced, although the ions eventually reached a temperature similar to that 
observed for smaller particle number. As statistical fiuctuations are suppressed only with 
V^, simulations with low density ratio and a sufficiently large number of particles per cell 
are prohibitively ex pensive, even in two dimensions. Other computational techniques like 



Sf-PlC simulations ( ISydoralll999l ) may be better suited for that purpose. 



Also, a substantial part of the initial heating and small-scale electric fields is related 
to the initial electrostatic Buneman-type instability, the effect of which depends on the 
drift velocity of the background electrons. Because the electrons carry the return current 
that compensates the cosmic-ray current, their drift velocity is Vd = VshNcR/N^. We must 
keep the shock speed Vsh (the cosmic-ray drift velocity) high to enable a sufficiently fast 
growth of the non- resonant streaming instability. We can vary the density ratio Ncn/Ne 
somewhat using particle splitting, but that also adversely affects both the growth rate and 
the wavelength of the non-resonant streaming instability. We have run three test simulations 
with Ncn/Ni = 1/10 and one with Ncr/Ni = 1/30 to further explore this issue. We find that 
the plasma temperature, or more precisely the average random kinetic energy per particle, 
is indeed significantly reduced by about an order of magnitude throughout the simulation. 
Nevertheless, the buildup of magnetic turbulence is unchanged compared with the earlier 
experiments. The plasma temperature is still unrealistically high, probably on account of 
heating in small-scale electric fields that arise from statistical charge-density fiuctuations. 



-17- 



4. Summary and Discussion 

4.1. Summary of Our Results 

The generation of magnetic-field turbulence by cosmic-ray ions drifting upstream of SNR 
shocks has been studied using two- and three-dimensional PIC simulations for a variety of 
parameters. Turbulent field is indeed generated in this process, but the growth of magnetic 
turbulence is slower than estimated in the literature, and the turbulence is of a different 
nature: we observe a modified filamentation of the ambient plasma, but not the cosmic rays, 
in contrast to the parallel wave found in quasi-linear calculations. The filamentation and 
formation of cavities in the background plasma was also observed in recent MHD simulations. 

The amplitude of the field perturbations saturates at approximately the amplitude of 
the homogeneous upstream field. The energy density in the turbulent field is also always 
much smaller than the plasma kinetic energy density. This suggests that the efficiency of 
magnetic-field generation through this mechanism may not be sufficient to account for the 
strong magnetic-field amplification invoked for some young SNRs and also leaves open the 
question whether or not diffusive particle acceleration at SNR shocks can produce particles 
with energies beyond the "knee" in the cosmic-ray spectrum. 

The backreaction of the magnetic turbulence on the particles leads to an alignment of 
the bulk ffow velocities of the cosmic rays and the background medium in all our simulations. 
This is precisely what is making up a cosmic-ray modified shock: the upstream fiow speed 
is continuously changed by the cosmic rays, so the compression ratio of the actual shock of 
the thermal plasma is moderate, whereas the overall compression ratio can be large. The 
new and surprising result of our simulations is that we accomplish this without significant 
magnetic-field amplification and with instability modes different from those invoked for the 
purpose in the recent literature. 



4.2. Comparison with Published MHD Simulations 



How can we understand the differences between the results of the analytical solutions, 
the MHD simulations, and our PIC study? First of all, not too much is different: the fila- 
mentation mode observed in our simulations is similar to the cylindrical filamentary mode 
that was analytically described in B0 5 and in its ev o lved, t hen isotropic form observed in the 
MHD simulations of B04, BOS, Z08, iReville et al.l (j2008al ). The growth rate, or expansion 
rate, of the filaments is somewhat smaller than expected for isolated filaments, probably be- 
cause neighboring filaments produce between them magnetic fields of opposite orientation. 
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thus cancelling part of the magnetic field that one would calculate for an isolated current 
filament. We do not observe a filamentation of the cosmic rays, probably because through- 
out the simulations the cosmic-ray Larmor radius, rcRg-, remains larger than the spatial 
scale of the magnetic turbulence. As we observe the magnetic-field growth to saturate and 
the cosmic-ray drift relative to the background plasma to disappear, there is no significant 
cosmic-ray current remaining that could initiate a further amplification of the magnetic field 
at later times, beyond the termination of our simulations. 

The MHD simulations of B04, BOS, and Z08 show the magnetic-field growth to an 
amplitude much higher than the initial homogeneous field. However, they assume the cosmic- 
ray current to be constant throughout the simulations. In our study the cosmic-ray current 
changes and is strongly reduced in the non-linear phase, when the magnetic-field growth 
saturates. Because in the MHD simulations the cosmic-ray current is constant in time and 
uniform in space, these simulations consequently cannot capture this backreaction. In the 
test run I* (see §3. 51 and Table 1), in which the backreaction of cosmic rays has been excluded 
by increasing cosmic-ray particle mass, we observe the peak in the average magnetic-field 
strength more than a factor 2 higher than with cosmic-ray backreaction (run I). In this test 
simulation the saturation of the field amplification arises from the bulk acceleration of the 
ambient plasma which is still permitted. In the MHD simulations, it is unclear to what 
extent momentum is transferred to the background plasma, but in any case that process is 
likely suppressed because the energy and momentum fiux of the cosmic rays is not accounted 
for. 

B04, BOS, and Z08 describe the magnetic-field structure in the MHD simulations at 
various stages of the evolution. The non-linear turbulence appears fairly isotropic and bears 
some resemblance to the structure shown in Figs. [5] and [TJ although the MHD results reveal 
no structures visibly extended in the drift direction. The MHD simulations show much 
more structure on the smallest scales, most of which are prominent at late stages when the 
turbulence in our PIC simulations has already saturated. The one-dimensional spectra of 
the perpendicular magnetic field in the MHD simulation of Z08 indicate strong growth on all 
spatial scales after the turbulent field becomes stronger than about 10% of the homogeneous 
field. From that time on the velocity fiuctuations in the background plasma are also of the 
same order as, if not stronger than, the adiabatic sound speed, indicating that weak shocks 
can be formed. In our simulations shocks would be resolved, which may explain why the 
MHD simulations show more structures on the smallest scales. 

The discussions in B04, BOS, and Z08 do not indicate that a parallel planar mode is 
indeed initially observed in the MHD simulations. Also, neither magnetic-field spectra in k± 
nor spectra of the parallel component of the magnetic field are presented. Our discussion 
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of the differences between the results of these MHD simulations and our PIC study must 
therefore remain limited in scope. We stress, however, that the saturation of the turbulence 
growth in our simulations arises from the backreaction of the magnetic field on the bulk 
velocities of cosmic rays and background plasma. The higher magnetic-field saturation level 
in the MHD simulations is most likely a direct result of the assumption of a constant cosmic- 
ray current. A strong magnetic-field amplification to amplitudes 5B 3> Bq has yet to be 
demonstrated. 



4.3. Comparison with Analytical Calculations 

The turbulence observed in our simulations reflects analytical results for filamentation 
modes, e.g., the derivation by BOS. A discrepancy only exists as far as the parallel plane- 
wave mode is concerned. All published analytical calculations agree that in quasi-linear 
treatment a parallel, purely growing mode should dominate. What might be the reasons 
why we do not see this mode in our PIC simulations? We have already determined that the 
plasma temperature is probably not responsible for this discrepancy. One possibility is that 
the mode exists initially, but is invisible in the noise, and changes its character quickly so 
that we observe the filamentation, which was already described in BOS. In fact, the parallel 
plane-wave mode and the filamentation of the ambient plasma are expected to show a similar 
growth rate. This interpretation would be in line with the magnetic-field structure observed 
in the MHD simulations before the turbulent magnetic field reaches the amplitude of the 
homogeneous field, Syo, which also does not clearly show a parallel mode. 

The most likely reason is that our choice of parameters may not reflect one of the as- 
sumptions made in the analytical treatments, nam ely that the frequen cy of the perturbations 



be much smaller than the ion gyrofrequency (e.g. iReville et al.l 120071 ). whereas we typically 
have to choose parameters for which the theoretically expected growth rate is similar to the 
ion gyrofrequency, i.e. ^ui ~ According to the quasi- linear results, 



'^u Vsh N^ 



CR ^^Q^ 



ni 2VANi ' 

implying that if either the shock speed or the cosmic-ray density is too high, the assumption 
of a small frequency is no longer justified, and the character of the instability may be different. 
However, the observed growth rate of the filamentation mode in our simulations is an upper 
limit to the growth rate of any parallel mode that we do not see, so in fact all our simulations 
have an observed growth rate '^Uobs ^ O.SQi and some simulations show '^Uobs < 0.2 fij. 

If So; ^ Qi must be strictly maintained, then what are the constraints on the parameters 
so the parallel plane-wave mode can play a role? Assuming efficient Bohm-type diffusion. 
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the upstream plasma is swept up by the shock on a timescale crcRg/Sv'^f^, so efficient growth 
requires the growth time be substantially smaller than this: 

If conditions allow a rapid growth of the instability on a timescale much shorter than the 
shock-capture timescale, meaning if the instability is astrophysically relevant, then the in- 
stability will also evolve in an environment that is essentially not changed by the inflow of 
fresh material, implying that our simulation setup using periodic boundary conditions on all 
sides is appropriate. 

For young SNRs the forward-shock speed is ^ 0.01, so that only a narrow range of 
interesting parameters exists, unless we consider very high-energy cosmic rays for which the 
instability-driving current is small. We can use Eq. [inito rewrite the second relation in Eq. [TT] 
using UcR ^ NcR a/ 7ci? " 1 "^i as the cosmic-ray energy density, and Ubuik = \Nimi v'^^ 
as the bulk kinetic energy density of the upstream plasma: 

UcR > 12— Ubuik ~ 0.1 Ubuik. (12) 

Vsh 

Here the second relation applies to SNRs with small upstream plasma density Ni < 0.1 cm~'^, 
such as SN 1006, RX J1713-3946, or Vela Junior, for which the Alfven speed is va — 30 km/s 
for a reasonable magnetic-field strength i?o — 5 fiG. It is unlikely that the energy density in 
cosmic rays is much larger than the bulk energy density of the upstream plasma, suggesting 
that this instability can operate for only a few e-folding times in those SNRs. Remnants 
that expand into a high-density environment may be better suited for a strong growth of the 
non-resonant mode, but in those cases one has to consider ion-neutral collisions which can 



reduce the growth of the instability (IReville et al.ll2008bl ) 



The relation Eq. [12] may be a more severe constraint than the requirement 



hma.rcR,--^--^^^^»l (13) 



for the instability to exist in the ffist place (compare Eq. [T]), and should be used in addition 
to it. Equation 12 also shows that cosmic rays of very high energy are not necessarily better 
triggers of magnetic-field amplification. The relevance of cosmic rays in a certain energy 
band primarily depends on the energy density they carry. 

Our simulations used cosmic rays with a moderate Lorentz factor jcr = 2. It may be 
interesting to speculate how our results might change, had we been able to use a substantially 
higher cosmic-ray Lorentz factor, e.g., 'Jcr = 1000. A spatial re-organization of the cosmic 
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rays was not observed in our simulations, and higher-energy particles are more difficult to 
concentrate in certain locations, so the homogeneity of the cosmic-ray distribution would 
most likely not change. Leaving the energy density in cosmic rays unchanged, we would 
expect a much slower growth of turbulence on much larger scales. The condition Quj -C fij 
should always be on account of the small number of cosmic rays required to carry their 
energy density (see Eq. [10]), and so the parallel mode should be initially excited and then 
give way to the filamentation and formation of cavities that was observed in the MHD 
simulations and our PIC studies. For the backreaction and saturation we therefore do not 
expect fundamental deviations from the behavior that we saw in our simulations. We studied 
systems with turbulence growth on a variety of scales relative to the plasma scale, and we 
never saw a systematic variation in the saturation mechanism and level, that might indicate 
a dependence on the wavelength of the turbulence. Even though we could not possibly study 
a system with, e.g., TeV-ish cosmic rays, we therefore feel there is no reason to assume that 
in such a situation the bulk acceleration of the plasma and the cosmic rays would proceed 
in a different way. The saturation level of the turbulent magnetic field would then also 
be similar to what we find. What remains unclear is the actual equilibrium amplitude in 
a realistic astrophysical scenario, in which the instabilities operate under a competition of 
saturation through non-linear backreactions with driving through the influx of fresh material. 
Since the saturation level observed in our simulations corresponds to a situation in which 
the backreactions are relatively fast, also compared with the initial growth, it appears likely 
that even under continuous inflow of new plasma the equilibrium amplitude of the turbulent 
magnetic fleld does not exceed the peak values seen in the PIC simulations, so that 6B /5 ~ 1 
remains the most probably result. 

Part of the simulations have been performed on Columbia at NASA Advanced Super- 
computing (NAS). Partial support by the National Center for Supercomputing Applications 
under PIIY070013N is acknowledged, where we used the Tungsten and Mercury systems. 
The work of J.N. was supported by MNiSW during 2005-2008 as research project 1 P03D 
003 29 and The Foundation for Polish Science through the HOMING program, which is sup- 
ported by a grant from Iceland, Liechtenstein and Norway through the EEA Financial Mech- 
anism. K.-I.N. is supported by AST-0506719, HST-AR-10966.01-A, NASA-NNG05GK73G, 
and NASA-NNX07AJ88G. 



A. Transformations of the Particle Distributions and their Momenta 

It is convenient to use cylindrical phase-space coordinates p\\, p±, and (p, where the 
differential volume element is d^p = dp\\ dp±p±d(j) = \dp\\ dp\ dcp, with the axial component 
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p\\ defined, in this case, by the direction of a Lorentz transformation. In our discussion the 
azimuthal angle can be ignored, so it is sufficient to consider gyrotropic distribution functions 



dn 



TT d?x dpw dp] 



(Al) 



If the distribution function is known in some inertial reference frame K, then in a second 
reference frame K' moving at velocity (3ce\\ with respect to K, the density, bulk velocity, 
and total energy density of the particles are given by 



n' ^ Jd'p'fip') 
V = l|dVv'/'(p) 
^' ^ I d'p'E'fip'), 



(A2) 
(A3) 
(A4) 



where primed quantities are measured in frame K'. We exploit the invariance of / and p± 
under the Lorentz transformation to rewrite Eq. IA2I as 



oo oo 



n' = 71 J dp\ j dpi f {pii ,p±) . 



-oo 



By the appropriate change of coordinates, the effect of the Lorentz transformation can be 
confined to the Jacobian of the coordinate transformation: 



oo oo 

„'^./*„/*i/(p,.PJ 



dp\\ 



(A5) 



px=const 



Now p'l = 7 — j3 -^Ipj , where E (p) = wm^c^ + + p^j Then the Jacobian is 



dp\\ 



dp\\ 



7 1- 



/3c 



(A6) 



Px=const. 

If we consider a particle distribution that is isotropic, homogeneous, and monoenergetic with 
scalar momentum pq and density n in frame K, 

n 



/(P) 



2lTpo 



(A7) 



we will find that n' = jn. The Jacobian must also be used when setting up the simulations, 
because in their flow frame cosmic rays are supposed to follow a distribution according 
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to Eq. IA71 which must be transformed into the simulation frame, which is the upstream 
rest-frame ahead of the SNR shock. 

We can compute the bulk velocity (Eq. IA3p in a similar manner; we express the particle 
velocity v' in terms of the integration variables, and only the parallel component f y will 
survive integration. Now = {v\\ — /5c)/(l — Pv\\/c), and v\\ = p\\ c^/E. Then for isotropic 
and monoenergetic particles (Eq. IA7I) we may write Eq. IA3I as 

oo oo 

'fn J J h — (jp\\c 

-oo 



dp\\ 



{At 

n I =const 



When Eq. lASl is solved with the appropriate substitutions (Eqs. IA6l and R7|l . the bulk 
velocity is found to be 

V' = -/5ce||. (A9) 

Thus the velocity for the Lorentz transformation to the rest frame K of the particle distri- 
bution is equal to the bulk, or drift, velocity of the particles in another frame K' , such as 
the simulation frame. The particle distributions are not necessarily isotropic in any frame 
of reference, but Eq. IA8I nevertheless allows us to properly calculate the drift velocity of 
a particle population and to transform the particle distribution into the flow rest-frame, so 
that the anisotropy properties can be investigated. We can also distinguish bulk and random 
momentum, and the transformed distribution function in the flow rest-frame gives the par- 
ticle distribution in random momentum that carries information on heating of the particles 
and any deviations from the initially prescribed Maxwellians for the background plasma and 
monoenergetic distributions for the cosmic rays. 

The total energy density w' is easy to calculate as well, but it is useful to separate it 
into rest-mass, random and bulk kinetic energy densities as w' = n' mc^ + U^.^^ + t/buik- -^^r 
simplicity we define them as U^^^ = (w'/t) — n'mc^ and f/^uik = (7 ~ 1)'W^V7) as a more 
detailed calculation would yield for an isotropic distribution. Fig. [1] shows the random and 
bulk kinetic energy densities (without primes) for the main three-dimensional experiment, 
run A. Note that the drift velocities, and hence 7, continuously change throughout the 
simulation. 
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Fig. 1. — For the main three-dimensional simulation (run A) using ion-electron mass ratio 
mi/nif. = 10, we show the temporal evolution of the average energy density in electromagnetic 
fields (left panel) and particles (right panel), both normalized to the initial total kinetic 
energy in the system, Uko- Time is in units of the inverse growth rate for the most 
unstable mode as predicted in the analytical theory (Eq. 2). B"^ and El_ indicate the 
magnetic and electric field energy densities in the components perpendicular to the cosmic- 
ray ion drift direction, i.e. {By + B'^)/{2fio) and correspondingly for the electric field. Bj, 
SBp and -E^ are defined analogously, where involves only the turbulent component. 
In the right panel, the total energy density of cosmic rays, ambient electrons, and ions 
is split into bulk and random components as outlined in Appendix |Al The total kinetic 
energy density in cosmic-ray ions, Ucr, and its components are scaled by a factor of 10. For 
comparison, the right panel also shows the time evolution of the volume-averaged magnetic 
energy density, Um = {Bl + B^ + 52)/(2/io). 
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Fig. 2. — The time evolution of the Fourier power spectrum of the perpendicular magnetic- 
field component By in run A. The spectra in wavelengths along the drift direction A^,. = 2Ti/k\\ 
are shown in the left panel (a), whereas spectra for the direction perpendicular to both the 
drift velocity and By are displayed in the right panel (b). Note that at small wavelengths 
(A < 5A) strong filtering reduces the Fourier power to very small values. The spectra are 
normalized to the peak spectral density. 
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Fig. 3. — The temporal evolution of the electron-density power spectrum for run A. The 
spectra are set up as described in Fig. [2] and normalized to the peak spectral density. 



- 29 - 




0.0 0.2 0.3 0.5 0.6 0.8 

Fig. 4. — Ambient electron-ion density contrast N^, — Ni normalized to the initial ion density 
for run A at a;/ A = 500 and t ~ 107~^^, as in Fig. [51 The average value for the normalized 
density contrast is 1/3 on account of the excess electrons that are needed to neutralize the 
charge carried by the cosmic-ray ions. There are no systematic large-scale deviations from 
the mean, indicating that the electron and ion density distributions are nearly cospatial. 
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Fig. 5. — The left panel displays the magnitude and direction, indicated by the arrows, of 
the perpendicular magnetic field component B± = (By + B'^Y^'^ in the plane perpendicular 
to the cosmic-ray drift direction at the grid position x/A = 500 and time t ~ ^O'j^l^. B± 
is normalized to the amplitude of the homogeneous field i?||o. The right panel shows the 
density of ambient ions, Ni, in units of the initial ion density at the same location and time. 
The electron distribution follows that of ambient ions (see Fig. H]). 



- 31 - 




0.1 0.5 0.9 1.3 1.7 2.2 




100 200 300 400 500 600 700 800 900 

^ x/A 



0.01 0.36 0.72 1.08 1.43 1.79 

Fig. 6. — The magnitude and direction of the perpendicular magnetic field component B± = 
{By + Biy/"^ (bottom panel) and the ambient ion density iVj (top panel) in the plane of 
the cosmic-ray ion drift direction at y/A = 150 and t ~ ^^Imax- ^ comparison with Fig. [5] 
illustrates the appearance of oblique filamentary structures. 
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Fig. 7. — Snapshots of the time evolution of the ion density and the perpendicular magnetic- 
field component structures in the plane perpendicular to the cosmic-ray ion drift direction 
at x/A = 500. 




10 20 30 40 

t7mox 



10-^ 



10" 



10" 



10" 



b 



0.1 U 



CR 



U; / 



10 



-lOOOOOOOOOOOCCiOOOOOOO 



Run A* 



20 

t7mo. 



30 40 



Fig. 8. — Temporal evolution of the total energy density in electromagnetic fields and par- 
ticles for run A* (see Fig. [T]). 
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Fig. 9. — Temporal evolution of the Fourier power spectrum of the magnetic field component 
for run A* (see Fig. [2]). 
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Fig. 10. — Temporal evolution of the total energy density in electromagnetic fields and 
particles for the two-dimensional simulation with a realistic mass ratio mi /me = 500 (run 
E; see Fig. HI). 
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Fig. 11. — Temporal evolution of the total energy density in electromagnetic fields and 
particles for the two-dimensional simulation assuming \max = 150 and mi/rrie = 10 (run F; 
see Fig. [1]). 
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Fig. 12. — Temporal evolution of the Fourier power spectrum of the magnetic field component 
Bz for run F (see Fig. [2]). 



-35- 




Table 1. Simulation Parameters and Results 
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Note. — Parameters and selected results of the simulation runs described in this paper. Listed are: the system size 
(a; x y X in units of the grid cell size A, the run duration in units of i'^ax-> ion-electron mass ratio rrii/me, 
the density ratio of ambient ions and cosmic rays, the ion skindepth A^j in units of A, the wavelength of the theoretically 
expected most unstable mode Xmax (in units of A, sec Eq. 2), the Alfvenic Mach number = Vsh/va of the shock, the 
plasma magnetization as given by the electron plasma to cyclotron frequency ratio, the maximum growth rate imax/^pe 
(see Eq. 2), the actual measured growth rate 7 in units of imaxi the obliqueness 6 = Z{k,Vs) in degrees of the actual 
dominant turbulence mode, and the maximum amplitude of the perpendicular magnetic-field perturbations SB^°'^ relative 



to the homogeneous magnetic field. 

^Cosmic-ray particles' mass assumed in this run is 2xl0^mi. 
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